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Abstract 

In this work, we highlight the correspondence between two descriptions of a system 
of ultracold bosons in a one-dimensional optical lattice potential: (1) the discrete 
nonlinear Schrodinger equation, a discrete mean-field theory, and (2) the Bose- 
Hubbard Hamiltonian, a discrete quantum-field theory. The former is recovered from 
the latter in the limit of a product of local coherent states. Using a truncated form of 
these mean-field states as initial conditions, we build quantum analogs to the dark 
soliton solutions of the discrete nonlinear Schrodinger equation and investigate their 
dynamical properties in the Bose-Hubbard Hamiltonian. We also discuss specifics of 
the numerical methods employed for both our mean-field and quantum calculations, 
where in the latter case we use the time-evolving block decimation algorithm due 
to Vidal. 
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1 Introduction 



Ultracold atoms in optical lattices have recently been at the center of excit- 
ing research by both experimental and theoretical physicists alike. In these 
systems, experimentalists are given an unprecedented amount of control over 
system parameters, and theorists are able to accurately describe the physics 
with relatively simple, clean models. However, a full quantum many-body 
treatment of the problem is computationally challenging due to the exponen- 
tial growth of the Hilbert space with the system size. There do exist numer- 
ous numerical methods, such as quantum Monte Carlo [1] and density matrix 
renormalization group |2], that can accurately calculate the ground state of 
the governing quantum Hamiltonian. On the other hand, for systems not ex- 
hibiting strong correlations, one can safely employ an appropriate mean-field 
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theory and describe the many-body system semiclassically, resulting in a much 
more tractable mathematical problem. Also, before the realization of ultracold 
atoms, dynamical properties of the many-body lattice problem saw little atten- 
tion; however, far-from-equilibrium quantum dynamics has recently become a 
hot topic of research [3|1] . Numerical work in this area has been spurred by the 
advent of a newly available algorithm for simulating the quantum dynamics 
of one- dimensional lattice Hamiltonians [5]. 

In this paper, we study dynamical properties of ultracold bosons on one- 
dimensional (ID) optical lattices. Our presentation is organized as follows. 
First, we elucidate the correspondence between the mean- field and quantum 
many-body theories in the case of the tight-binding approximation on a lat- 
tice. Second, we discuss our numerical methods, including an implementation 
of the time-evolving block decimation algorithm to simulate the quantum dy- 
namics [5] . Finally, we simulate and analyze the time evolution of dark solitons 
under both the semiclassical and quantum equations of motion and highlight 
the observed differences between the two theories. 



2 Quantum-Mean Field Correspondence 

2.1 Mean-Field Theory 

The statics and dynamics of a weakly interacting Bose gas at zero temperature 
in free space, or in the more experimentally relevant geometry of a harmonic 
trap, are well- described by the Gross-Pitaevskii equation (GP), a.k.a., the 
nonlinear Schrodinger equation (NLS) [B]. For these simple geometries, even 
in the quasi-lD regime, such a mean-field approach is appropriate given that 
there is negligible depletion out of the condensed mode. However, quantum 
fluctuations cannot be ignored throughout time evolution for certain excited 
condensate states such as the dark soliton [7f5] . Mean-field theory treatments 
have been generalized to lattice geometries [H] in which case the lattice soli- 
ton solutions of the NLS have been mapped out in detail [TO|ITT] . However, 
such analysis assumes that the lattice height is sufficiently low so that the 
gas still exhibits nearly ideal Bose-condensation, i.e., negligible depletion out 
of the single boson configuration as assumed by the NLS. When using an 
unperturbed mean-field theory to describe the system, one ideally uses the 
continuous NLS with an external lattice potential; however, coupled-mode 
theory can be employed for a shallow lattice [11] or, in the other extreme, a 
single-band tight-binding approximation can be used for a deep lattice [T2] . 

In the latter case of the tight-binding approximation, assuming ID from here 
on, the full condensate wave function t) is expanded in a basis of localized 
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condensate wave functions (j){x — Xi) each centered at site i in the lowest Bloch 
band of the lattice. That is, t) = J2i 'ipi(t)(f){x—Xi). It is then assumed that 
the local condensate wave functions are sufficiently localized in each well, in 
which case we obtain the well-known discrete nonlinear Schrodinger equation 
(DNLS): 

ihdttpk = - J {'4'k+i + ^fc-i) + Ulipk^ipk + CfcV'fc, (1) 

where Vfc = ^fc(^) is the dimensionless c- number that weights the localized 
condensate wave function at the kih. lattice site, and is a local external po- 
tential different from the lattice potential. The coupling parameter J and the 
nonlinearity U can be written explicitly in terms of overlap integrals involving 
the system parameters and the localized condensate wave functions (f)[x — Xi) 



2.2 Quantum Many-Body Theory 



We now turn to the full quantum many-body description of the problem in 
the regime of optical lattice depths where the above mean-field tight-binding 
description is appropriate, and beyond. As was demonstrated by Jaksch et al. 
[T3] , a system of weakly interacting ultracold bosons loaded in an optical lattice 
potential is an almost perfect realization of the Bose-Hubbard Hamiltonian 
(BHH), a model introduced to the condensed matter community almost ten 
years earlier by Fisher et al. 



To derive the BHH, we start with the ID continuous many-body Hamiltonian 
in second quantization for two-body interactions: 



H = j dxi>\x) 



Kxt(a;j 



2m dx"^ 



^!{x) 



+ ^JdxJ dx'¥{x)'^^x')Vir,t{x - . (2) 

We then expand the bosonic field operator which destroys a particle at 

position X, in a lowest Bloch band Wannier basis as ^'(x) = J2ibiw{x — x,), 
where the operator bi is defined to destroy a particle in the localized Wannier 
wave function w{x — Xi). This step is analogous to expansion of the full con- 
densate wave function in a localized basis when discretizing the continuous 
NLS to obtain the DNLS. Then, as is done in the derivation of the continuous 
NLS, we assume the two-body interaction potential to be of the contact form, 
i.e., Viai{x — Xi) = g6{x — Xi), where g is proportional to the s-wave scattering 
length of the atoms. We also invoke the tight-binding approximation by as- 
suming that the lattice is deep enough to obtain sufficiently localized Wannier 
functions. This allows us to discard all terms except those involving nearest- 
neighbor hopping and on-site interactions. After making these assumptions. 
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we arrive at the familiar BHH: 

M-l ^ ^ U ^ ^ 

H = -JJ2 (bl+ibi + h.c.) + - 51 - i) + E (3) 

i=l ^ i=l 1=1 

where 6j and 6| are destruction and creation operators at site i that obey the 
usual bosonic commutation relations, and nj = is the number operator 
which counts the number of bosons at site i. Equation assumes box bound- 
ary conditions on a lattice containing M sites. The coefficients J, U, and can 
be calculated exactly in terms of the localized single-particle wave functions 
and other parameters [13]. J is the nearest-neighbor hopping coefficient, U is 
the on-site interaction energy, and is an external potential. The ratio U/J is 
an important parameter that determines the relative contribution from each 
term. For a shallow lattice, the hopping term dominates, whereas for a deep 
lattice the interaction term dominates. However, we note that U is not com- 
pletely dependent on the lattice geometry since the s-wave scattering length 
can be varied independently via a Feshbach resonance. 

In deriving the BHH, one makes very similar assumptions to those made when 
discretizing the continuous NLS on a lattice to obtain the DNLS. Specifi- 
cally, both derivations invoke a lowest Bloch band tight-binding approxima- 
tion. However, in the latter case, a single configuration of bosons is assumed 
from the onset. The full quantum treatment allows for quantum depletion 
out of the condensate mode and thus can describe the system in the strongly 
interacting regime. 

A general pure state of the full many-body quantum system can be written 
as a complex linear superposition of states, each with a well-defined number 
of particles in each Wannier state: 

d-l 

I*) = H Cn,n2-nM\nin2 ■ ■ ■ Um) , (4) 

ni,n2,...,nM=0 

where is the number of particles at site k. For obvious computational rea- 
sons, we truncate the local Hilbert space at local dimension d, i.e., we restrict 
the occupation of each Wannier state to contain at most d — 1 bosons. The 
Hilbert space containing all pure states of the full many-body system is thus 
of dimension d^'^ which becomes prohibitively large for large systems. For ex- 
ample, even for d = 2 we can only simulate M = 12 or 13 lattice sites on a 
single PC without further refining the numerical algorithm. We overcome this 
difficulty with use of the time-evolving block decimation routine which will be 
discussed in Sec 13.11 
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2.3 Discrete Mean-Field Theory From Discrete Quantum Many-Body Theory 



Next, we show how the DNLS can be recovered from the BHH. The destruction 
operator at site k can be evolved in time in the Heisenberg picture according 
to ihdtbk = [bk,H]- After computing the commutators, we arrive at 

ihdtbk = - J {bk+i + bk-i) + Ubkb\bk + ekh- (5) 

We can then take the expectation value of Eq. ^ to obtain an equation of 
motion for the order parameter (6^). The DNLS is recovered exactly if the 
expectation value is taken with respect to a product of atom-number Glauber 
coherent states. That is, for full many-body states of the form 

I*) = l^fc)' where \zk) = J2 ^1^)' (6) 

k=i n=o V ri! 

we obtain the DNLS for the equation of motion governing the coherent state 
amplitude Zk = {bk). 

ihdtZk = -J {zk+i + Zk-i) + U\zk\^Zk + ekZk. (7) 

As discussed in Sec. 12.21 for numerical calculations we must truncate the local 
Hilbert space to a finite dimension d, in which case the on-site coherent states 
of Eq. ([6]) become truncated coherent states: 

M 1^ |2 d-l 

1^) = l^fc)' where \zk) = Xde'^ ^I'^)' (8) 

k=l n=0 V J^' 

and A/rf is a normalization factor. 

The coherent states of Eq. ([6]) are known to well-describe the ground state of 
the BHH for J ^ U in the limit of an infinite number of sites M and particles 
at fixed filling N/M [13]. It is in this regime that quantum depletion can be 
safely neglected and Eq. ([7]) is an accurate description of the system. However, 
the lattice must still be deep enough so that the single-band tight-binding 
approximation is still valid. In Sec. HI we use the truncated coherent states of 
Eq. (IH]) to create nonequilibrium initial quantum states in the BHH that are 
analogs to the dark soliton solutions of the DNLS. 
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3 Numerical Methods 



3.1 Time- Evolving Block Decimation Algorithm 



The time-evolving block decimation (TEBD) algorithm was first introduced 
in 2003-2004 by Vidal [T6|f5] in the context of quantum computation. Soon 
thereafter, Daley et al. [17] and White and Feiguin [18] translated the algo- 
rithm into more the familiar density matrix renormalization group (DMRG) 
language and showed that TEBD is equivalent to a time- adaptive DMRG 
routine. Here, we summarize our implementation of TEBD as applied to the 
BHH. 



3.1.1 The Vidal Decomposition 



The Vidal prescription is to first rewrite the coefficients in Eq. (j4]) as a product 
of M tensors {EM} and M - 1 vectors {A^}: 



E 



F[l]«i\[l]F[2]n2 \[2]-p[3]n3 . . ,-p[M]nM (n\ 
oi,...,aA/_i=l 



There does exist a procedure for determination of the Fs and As given known 
coefficients of an arbitrary state; however, this is not generally useful because 
one does not typically have access to each component of the c/^^-dimensional 
vector This procedure would require a Schmidt decomposition (SD) at 
every bipartite splitting of the lattice, where x ^ c/L*^/^-! is the number of 
Schmidt vectors retained at each splitting. The Schmidt number xs^ i-e-? the 
number of Schmidt basis sets required for an exact representation of the state 
at each cut, is naturally a measure of global entanglement between the lat- 
tice sites [T6IT9] . The decomposition iQ is thus appropriate when |\E') is only 
slightly entangled according to the Schmidt number, in which case it is com- 
putationally feasible to take x ~ Xs- 



3.1.2 Two- Site Operation 

One of the reasons why this decomposition is useful is that it allows for efficient 
application of two-site unitary operations. Let us consider a two-site unitary 
operation V = acting on sites £ and £+1. First, 

we write |\E') in terms of Schmidt vectors for the subsystems [1 ■ ■ — 1] and 
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[£ + 2---M]: 



by invoking Eqs. (13) and (14) of [TB], where 

on^n^+i ^ \[^-l]T-Mn, , [f] -pff +l]„,_^i , [£+1] /..n 
oil 

and G {1,2,..., x}. Note that this definition of the tensor differs from 
an analogous construct in [16] which is also denoted G in that work. We are 
up to this point assuming that we know the decomposition ^ of |\Ef), and 
hence we also know all elements of 6. However, by writing |\E') in the form of 
Eq. ffTOj) we can easily write the updated state after the application of V as 

where can be written in terms of the updated tensors r^^' and rt^+^l and the 
updated vector A'^^: 

In practice, a given two-site operation is performed as follows: (1) form B 
from current Fs and As [Eq. ( ITTi) ]: (2) update 6 by applying V to obtain 
6 [Eq. (fT3!) ]: (3) reshape 6 from a 4-tensor to a (xc?) x (xc?) matrix; (4) 
perform a singular value decomposition (SVD) on this matrix retaining only 
the largest x singular values A|^|; and (5) divide out the previous values of 
A^^"^' and A^^''"^] in order to compute F^^' and F'^"^^' from the matrix elements 
obtained via the SVD. The most expensive computational steps are (1), the 
formation of 0, and (2), the update of after the apphcation of V. The 
former requires 0{d?y^) elementary operations, whereas the latter requires 
0{d^'X^) elementary operations; hence, our overall two-site operation scales as 
C[max(d2;^3,dV)]- 



3.1.3 Real Time Evolution 

The BHH is a sum of one- and two-site operations, but the terms multiplying J 
in Eq. ([3]) do not all commute, so the time evolution operator g-*^*/'* does not 
directly factor into a product of one- and two-site unitary operations. However, 
because the BHH only links nearest neighbors, we write H = Hodd + -f^even, 
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where 



odd 



^ E(^l+i^^ + h.c.)+ E 



odd 



odd 



J E (^lA + h.c.)+ E 



u 



-hi{hi - 1) + eifii 





and 



(14) 
(15) 



Each term within both -ffodd and i^cvon commute even though [i/odd, -f^cven] 7^ 0. 
It is then convenient to utihze a Suzuki- Trotter approximation of the time evo- 
lution operator for small time steps 5t. Specifically, we employ the second-order 
expansion: e"*^*^*/^ Q-^iioddSt/2h^-iH^^cnSt/h^-tHaa^5t/2h^ where each exponen- 
tial factor can be factored into a product of two-site unitaries. Even though 
the terms involving n are one-site operations, we still treat them as two-site 
operations by appropriate tensor products with the identity operator. In prac- 
tice, we build d^-dimensional matrix representations of H for each lattice link 
and diagonalize these matrices to obtain matrix representations of the two-site 
unitary operators. Then, in conjunction with the Suzuki- Trotter expansion, we 
employ the two-site operation procedure outlined in Sec. 13.1.21 on an initial 
decomposed configuration |\E') 0{M) times for each of tf/6t total time steps, 
updating the decomposition at each step. It is straightforward to calculate 
single-site observables, e.g., the expectation value of the number operator (rifc), 
and two-site observables, e.g., the one-body density matrix (h\hj), by using the 
partial trace to calculate the reduced density matrix of the subsystem of inter- 
est. For example, to calculate (blbj), we first compute pij = trfc^jj|\Ef)(\E'| using 
Eq. Q for |\E') and then use (blbj) = tT{b^^b pij). Overall, our implementation 
of the TEBD algorithm scales as 0[M^-^ max((i^x^5 d'^X^)]- 



3.1.4 Sources of Error and Convergence Properties 

The TEBD algorithm makes two important approximations: (1) the retention 
of only the x most heavily weighted basis sets during a given two-site oper- 
ation (see Sec. 13.1.21) . and (2) the Suzuki- Trotter representation of the time 
evolution operator (see Sec. 13.1.31) . For the latter case, we find that for the 
results presented in Sec. |l]it is sufficient to use time steps of size 6t = 0.01 h/ J 
to obtain converged results. The former approximation is more subtle as its 
accuracy is directly related to the amount of entanglement present in the sys- 
tem. Specifically, in Sec. HI we time evolve mean-field initial states [see Eq. ([8])] 
for which x = 1 is sufficient for exact representation; however, unitary time 
evolution increases entanglement between sites. To ensure that our choice of x 
is sufficient, we run equivalent simulations with increasing values of x and look 
for convergence of calculated observables, e.g., average local number (fik). It 
is important to point out that, owing to the local nature of the expansion (Q, 
accurate calculation of nonlocal observables, e.g., off-diagonal elements of the 
single-particle density matrix, converge more slowly with respect to x- For the 
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observables and time scales presented in Sec. HI our results are converged for 
the specified values of x = 45, 50. We also note that the fidelity of truncation 
at X eigenvalues can be quantified by the sum of the non-retained eigenvalues 
after application of a two-site operation. This quantity can be interpreted as a 
measure of the amount of entanglement generated by the two-site operation. 
Typically, for a single two-site unitary, we find this truncation error to be less 
than or on the order of 10~^. 



3.2 Constrained Imaginary Time Relaxation in DNLS 



3.2.1 Fundamental Dark Soliton Solutions 

In Sec. m we use the TEBD routine to simulate the quantum evolution of the 
dark soliton solutions of the DNLS by using truncated coherent states as initial 
configurations. This requires knowledge of the set of coherent state amplitudes 
{zk\ corresponding to a discrete dark soliton. Using a standard Crank-Nicolson 
scheme for the time-stepping procedure, we calculate the standing dark soliton 
solution of the DNLS by performing constrained imaginary time relaxation on 
Eq. ([7]) with efc = 0. Specifically, we take an initial condition of form Zk = mxk 
where Xk is the position of the fcth site and x = is the center of the lattice 
and normalize the solution to A^'dnls = I^fcli kfcP Qcich. step of imaginary 
time. The stationarity of the solution is tested by subsequent evolution in real 
time. 



3.2.2 Density and Phase Engineering of Gray Solitons 



We also consider the case of two solitons moving toward one another at fi- 
nite velocity. These initial conditions are obtained via the methods of density 
and phase engineering for soliton creation [20] as applied to the DNLS. We 
first perform imaginary time relaxation on a uniform initial condition with an 
external potential of the form 



ek = Vo{ exp 



2a1 



-|- exp 



2a2 



(16) 



to dig two density notches each centered at distance ^ from the center of the 
lattice. Next, we imprint an instantaneous phase of the form 



dk = M 



\ — tanh 




H — tanh 


"2(xfc-0' 


-I 


I 2 


^9 


2 







(17) 



which gives the solitons equal-and-opposite initial velocities toward the center 
of the lattice. Phonon generation is minimized by appropriately tuning the 
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width (j0 of the phase profiles to the soliton depth as determined by Vq in the 
density engineering stage. 



4 Time Evolution of Quantum Solitons 

With an initial DNLS configuration {2:^(0)} obtained either via the proce- 
dure outlined in Sec. 13.2.11 for a single standing soliton or the procedure 
outhned in Sec. 13.2.21 for two colliding sohtons, we then build a product 
of truncated coherent states according to Eq. ([8]) for input into the TEBD 
quantum simulation routine. The Vidal decomposition ([9]) of a product state 
1^) = {g)^li (E^;io dlVfc)) is trivial to compute: 

where for the case of truncated coherent states cl'^^ = AO e~'^*'^/^4^- 

An extensive discussion of resuhs obtained using the above methodology to 
create nonequilibrium dark sohton initial states in the BHH is presented in 
Refs. pT|22] . We will summarize those results here. Being direct analogs of 
mean-field solitons, the initial conditions analyzed below do not conserve to- 
tal particle number, although quantum evolution does conserve total average 
particle number. That is, in this section, we consider the quantum many-body 
evolution of mean-field- like solitons and refer to these structures as quan- 
tum solitons. However, we stress that neither the discrete mean-field theory 
(DNLS) nor the corresponding quantum theory (BHH) are integrable systems, 
so these are not solitons in the mathematically rigorous sense. It is possible to 
density and phase engineer dark soliton states directly in the BHH that are 
eigenfunctions of the total number operator. In Ref. [22], we use a number- 
conserving version of the TEBD routine to generate and analyze the quantum 
dynamics of such states. Although some observables behave differently in this 
case, the conclusions reached are generally the same. 

4.1 Standing Solitons 

The DNLS assumes a single configuration of bosons, i.e., bosons are only al- 
lowed to occupy one single-particle orbital. However, for an M mode system, 
a full quantum treatment will permit bosons to occupy any of the M permissi- 
ble modes. For the system sizes accessible to explore with the TEBD routine, 
a standing soliton initial DNLS configuration exhibits a finite lifetime due to 
quantum effects indescribable by mean-field theory. Most notably, quantum 
evolution causes significant quantum depletion out of the initial dark soliton 
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Site Index k 



Fig. 1. Density measures for a standing quantum soliton. Quantum evolution of (a) 
average particle number, (b) condensate wave function [23], and (c) order parameter 
versus position and time for a standing DNLS dark soliton initial configuration. The 
dashed line in (a) indicates the 1/e decay time of the order parameter norm Ni^. 



configuration into higher order orbitals which fill in the soliton density notch, 
where as discussed in [21], the natural orbitals of the system are defined as 
the eigenf unctions of the one-body density matrix (h\hj) [23] . We find that the 
soliton lifetime is closely correlated to the the growth in quantum effects such 
as a decay in the order parameter norm Nb = Efcli |(&fc)P and growth in the 
generahzed entropy Q = ^ ~ jj X]fc=i ^^iPl) ■ Shown in Fig. [Tjis the evo- 
lution of density measures of a quantum soliton with parameters uU/ J = 0.35, 
z/ = 1, M = 21, d = 7, x = 45, where u = N^jsils/M is approximately equal 
to the average filling. 



4-2 Soliton- Soliton Collisions 



In Fig. [21 we display the quantum evolution of two colliding dark solitons 
where the initial conditions were obtained by density and phase engineering 
in the DNLS as summarized in Sec. I3.2.2[ Here, if the decoherence time, as 
measured by the decay time of the order parameter norm, occurs before or near 
the collision time, then there is a loss in elasticity of the soliton collision. For a 
fixed value of the effective nonlinearity uU / J, we can independently tune the 
decoherence time by changing the filling u without altering the initial density- 
phase profile of the sohtons. In Fig. [2], we depict this effect in three separate 
simulations with parameters uU/J = 0.35, M = 31, x = 50, Vq/J = 0.4, 
a^/a = 1, A6 = O.Sir, ag/a = 2, ^/a = 6 at filling factors u = 1,0.5,0.1, where 
a is the lattice constant. 
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Site Index k Site Index k Site Index k 



Fig. 2. Quantum soliton collisions and decoherence-induced inelasticity. Average 
particle number for two colliding quantum solitons at filling factors (a) u = 1, (b) 
u = 0.5, and (c) u = 0.1. The collision elasticity is decreased when the decoherence 
time (dashed lines) occurs at or before the time of collision. 



5 Conclusion 

In conclusion, we have summarized the derivations of the DNLS both from 
the continuous mean field and from the continuous quantum field. In both 
cases, we invoke a single-band tight-binding approximation; however, the lat- 
ter derivation from the quantum many-body perspective is more insightful 
as it is based on single-particle physics with single-orbital occupation not as- 
sumed until the end. Using Vidal's TEBD routine to propagate dark soliton 
DNLS configurations forward in time according to the BHH, we have shown 
that quantum fluctuations give dark solitons a finite lifetime and induce an 
inelasticity in soliton-soliton collisions. For a more extensive analysis of these 
results, we refer the reader to Refs. [2T1I22] . 

We thank Charles Clark, Ippei Danshita, and Jamie Williams for useful dis- 
cussions. This material is based upon work supported by the National Science 
Foundation under Grant No. PHY-0547845, as part of the NSF CAREER 
program. 
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